Pangolin ID: P09
Periodicity analysis
Peter Kamerman
29 Mar 2019
This script includes analyses on whether there is evidence of circadian variation in body temperature in pangolins. The analysis is for one animal (P09) and includes chronograms, actograms, periodgrams, and autocorrelation.
Import and summarise data
Import
# Body temperature data
data <- read_rds('data-cleaned/pangolin/P09-cleaned.rds')
# Globe data
globe <- read_rds(path = 'data-cleaned/pangolin/globe-temp.rds')
# Import photoperiod data
photo <- read_rds(path = 'data-cleaned/pangolin/photo-period.rds')Summarise
kable(skim(data))Skim summary statistics
n obs: 98783
n variables: 6
Variable type: character
| variable | missing | complete | n | min | max | empty | n_unique |
|---|---|---|---|---|---|---|---|
| date | 0 | 98783 | 98783 | 10 | 10 | 0 | 343 |
| ID | 0 | 98783 | 98783 | 3 | 3 | 0 | 1 |
| period | 0 | 98783 | 98783 | 13 | 13 | 0 | 4 |
| time | 0 | 98783 | 98783 | 8 | 8 | 0 | 288 |
Variable type: numeric
| variable | missing | complete | n | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|---|
| body_temp | 0 | 98783 | 98783 | 34.83 | 0.49 | 31.12 | 34.5 | 34.74 | 35.14 | 38.34 | ▁▁▁▇▇▁▁▁ |
Variable type: POSIXct
| variable | missing | complete | n | min | max | median | n_unique |
|---|---|---|---|---|---|---|---|
| date_time | 0 | 98783 | 98783 | 2017-04-04 | 2018-03-12 | 2017-09-22 | 98783 |
Process data
#-- Generate separate day, month, and year columns --#
data %<>%
mutate(day = day(date),
month = month(date,
label = TRUE),
year = as.character(year(date)))
#-- Smooth body temp data --#
## Rolling smooth over 60 minutes (alignment = centered)
data %<>%
mutate(temp_smooth = round(rollapply(data = body_temp,
FUN = mean,
width = 12,
partial = TRUE,
align = 'center'), 2))
#-- Define animal ID --#
animal <- data$ID[[1]]
#-- Photo-period: generate separate day, month, year columns --#
photo %<>%
mutate(day = day(date),
month = month(date,
label = TRUE,
abbr = FALSE),
year = as.character(year(date)))
#-- Photo-period: spread time of sunrise/sunset
photo_wide <- photo %>%
spread(key = event,
value = time) %>%
group_by(date) %>%
# Shift sunset up one row
mutate(sunset = lead(sunset)) %>%
ungroup() %>%
# Remove empty rows
filter(!is.na(sunset)) %>%
arrange(date)
#-- Photo-period: filter to length of 'data' --#
photo_wide %<>%
filter(date %in% data$date)Daily summary line plots
#-- Body temp --#
data_summary <- data %>%
# Group and summarise
group_by(date) %>%
summarise(mean = mean(body_temp, na.rm = TRUE),
min = min(body_temp, na.rm = TRUE),
max = max(body_temp, na.rm = TRUE)) %>%
ungroup() %>%
# Add dummy time column and rejoin with data
mutate(time = '12:00:00') %>%
left_join(data) %>%
# Clean-up
select(date, period, mean, min, max) %>%
# Add body temp label
mutate(type = 'Body temperature')
#-- Globe --#
globe_summary <- globe %>%
# Group and summarise
group_by(date) %>%
summarise(mean = mean(globe_temp, na.rm = TRUE),
min = min(globe_temp, na.rm = TRUE),
max = max(globe_temp, na.rm = TRUE)) %>%
ungroup() %>%
# Add dummy time column and rejoin with data
mutate(time = '12:00:00') %>%
left_join(globe) %>%
# Clean-up
select(date, period, mean, min, max) %>%
mutate(mean = ifelse(is.nan(mean),
yes = NA,
no = mean),
min = ifelse(is.infinite(min),
yes = NA,
no = min),
max = ifelse(is.infinite(max),
yes = NA,
no = max))
#-- Globe filtered --#
# Create list of globe temps limited to body temp range of the animal
globe_filtered <- data_summary %>%
select(date) %>%
left_join(globe_summary) %>%
# Add a globe temp label
mutate(type = 'Globe temperature')
# Bind the globe and body temp data together
plot_data <- data_summary %>%
# Bind
bind_rows(globe_filtered) %>%
# Convert date to ymd format
mutate(date = ymd(date))
#-- Plot --#
ggplot(data = plot_data) +
aes(x = date) +
geom_ribbon(aes(ymax = max, ymin = min), fill = 'wheat2') +
geom_line(aes(y = mean)) +
geom_line(aes(y = max), colour = 'wheat3') +
geom_line(aes(y = min), colour = 'wheat3') +
labs(title = str_glue('{animal}: Daily mean and range in body and black globe temperature'),
subtitle = 'Black line: daily mean \nShaded area: daily range',
x = 'Date',
y = expression(Temperature~(degree*C))) +
scale_x_date(date_breaks = '2 months',
date_labels = '%b %Y') +
facet_wrap(~ type,
ncol = 1,
scales = 'free_y') +
theme_bw(base_size = 12) +
theme(axis.text.x = element_text(angle = 45,
hjust = 1),
panel.grid.minor = element_blank())
Chronograms
As a preliminary assessment of whether there may be circadian osscilation in body temperature, we plotted chronograms of raw and smoothed (rolling 60-minute average) body temperature data for each month in a calendar year.
Year facetted by month
#-- Nest data by year --#
chronograms <- data %>%
group_by(year) %>%
nest()
#-- Plot chronograms by month within each year --#
chronograms %<>%
# Add dummy date and date-time for plotting
mutate(data = map(.x = data,
~ .x %>%
mutate(dummy_date_time = ymd_hms(paste0('2000-01-', day, ' ', time)),
dummy_date = ymd_hms(paste0('2000-01-', day, ' 00:00:00'))))) %>%
# Plot
mutate(chronograms = pmap(.l = list(data, year, animal),
~ ggplot(data = ..1) +
aes(x = dummy_date_time,
y = body_temp) +
geom_vline(aes(xintercept = dummy_date),
colour = '#2678B2') +
geom_point(shape = 21,
stroke = 0.2,
fill = '#FFFFFF') +
geom_line(aes(y = temp_smooth),
colour = '#FD7F28') +
labs(title = str_glue('{..3}: {..2} monthly chronogram'),
subtitle = 'Open circles: 5-minute temperature data \nOrange line: Smoothed temperature data (60-minute rolling mean, center aligned) \nBlue vertical lines: Midnight of each day',
x = 'Day of month',
y = expression(Body~temperature~(degree~C))) +
scale_x_datetime(date_breaks = '2 days',
expand = c(0, 0),
date_labels = '%d') +
scale_y_continuous(limits = c(30, 38)) +
facet_grid(month ~ .,
drop = FALSE) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank())))
# Print plots
walk(chronograms$chronograms,
~print(.x))
Actograms
We then refined the visual inspection of the data for evidence of a circadian rhythm of body temperature by constructing monthly actograms that plotted difference in body temperature (averaged across 48 30-minute bins per day) from the daily mean body temperature for each day of a given month.
# Prepare photo-period data for plotting
photo_actogram <- photo_wide %>%
mutate(sunset = ymd_hms(str_glue('2000-01-01 {sunset}')), # Use a dummy date
sunrise = ymd_hms(str_glue('2000-01-01 {sunrise}'))) %>%
select(year, month, day, sunrise, sunset) %>%
# Convert day to an ordered factor
mutate(day = factor(day,
levels = 1:31,
ordered = TRUE)) %>%
mutate(month = as.character(month)) %>%
group_by(year, month) %>%
nest() %>%
rename(photo_period = data)
# Added later to fix cbind error
extra <- tibble(year = c('2018', '2018'),
month = c('February', 'March'),
photo_period = list(tibble(day = NA,
sunrise = NA,
sunset = NA),
tibble(day = NA,
sunrise = NA,
sunset = NA)))
photo_actogram %<>%
bind_rows(extra) %>%
mutate(photo_period = map(.x = photo_period,
~ .x %>%
mutate(day = factor(day,
levels = 1:31,
ordered = TRUE),
sunrise = ymd_hms(sunrise),
sunset = ymd_hms(sunset))))
#-- Nest data by year --#
actograms <- data %>%
mutate(month = month(date_time,
label = TRUE,
abbr = FALSE)) %>%
mutate(month = as.character(month)) %>%
group_by(year, month) %>%
nest()
#-- Join actogram and photo_actogram --#
actograms %<>%
bind_cols(photo_actogram[ , 3])
#-- Plot chronograms by month within each year --#
actograms %<>%
mutate(data_binned = map(.x = data,
~ .x %>%
# Create 30-minute bins
## Extract minutes
mutate(time_min = minute(date_time),
time_hour = hour(date_time)) %>%
## Re-level minutes to 30 minute intervals
mutate(time_30 = ifelse(time_min > 0 & time_min < 35,
yes = '00:00',
no = '30:00'),
time_30 = str_glue('{time_hour}:{time_30}')) %>%
# Bin by 30-min interval and calculate
# average temperature per bin
group_by(day, time_30) %>%
summarise(temp_binned = mean(body_temp, na.rm = TRUE)) %>%
# Calculate delta bin and mean daily temperature
group_by(day) %>%
mutate(temp_delta = temp_binned - mean(temp_binned)) %>%
ungroup() %>%
# Make time_30 a date-time object
mutate(time_30 = str_glue('2000-01-01 {time_30}'),
time_30 = ymd_hms(time_30)) %>%
# Generate a colour aesthestic for the
# direction of temp_delta
mutate(col_ours = ifelse(temp_delta < 0,
yes = 'Negative',
no = 'Positive')) %>%
# Convert day to an ordered factor
mutate(day = factor(day,
levels = 1:31,
ordered = TRUE)))) %>%
# Plot data
mutate(actograms = pmap(.l = list(data_binned, month, year, animal, photo_period),
~ ggplot() +
geom_rect(data = ..5,
aes(xmax = sunset,
xmin = sunrise),
ymax = Inf,
ymin = -Inf,
fill = '#FFFFFF') +
geom_col(data = ..1,
aes(x = time_30,
y = temp_delta,
fill = col_ours,
colour = col_ours)) +
labs(title = str_glue('{..4}: {..2} {..3} daily actograms'),
subtitle = 'Recordings averaged across 30-minute bins | Unshaded area: Sunrise to sunset',
x = 'Time of day',
y = expression(Delta~body~temperature~from~daily~average~(degree~C))) +
scale_x_datetime(date_breaks = '180 mins',
date_labels = '%H:%M',
expand = c(0, 0)) +
scale_y_continuous(breaks = 0) +
scale_fill_manual(name = expression(Direction~of~Delta),
values = c('#FD7F28', '#2678B2')) +
scale_colour_manual(name = expression(Direction~of~Delta),
values = c('#FD7F28', '#2678B2')) +
facet_grid(day ~ .,
drop = FALSE) +
theme_bw(base_size = 12) +
theme(axis.text.y = element_blank(),
panel.grid = element_blank(),
panel.background = element_rect(fill = '#CCCCCC'),
panel.spacing.y = unit(0, 'lines'))))
# Print plots
walk(actograms$actograms,
~print(.x))
Summary
The figure below is the monthly average of 5-minute difference in body temperature from mean daily temperature recordings for each year.
Please note that the sequence of months is September 2017 (orange) through to March 2018 (blue).
# Prepare photo-period data for plotting
photo_summary <- photo_wide %>%
# Extract first and last day of each month
group_by(month) %>%
filter(day == min(day) | day == max(day)) %>%
mutate(day = ifelse(day == max(day),
yes = 'D2',
no = 'D1')) %>%
ungroup() %>%
select(day, month, sunrise, sunset) %>%
unique(.)
photo_D1 <- filter(photo_summary, day == 'D1') %>%
mutate(sunrise = hms(sunrise),
sunset = hms(sunset))
photo_D2 <- filter(photo_summary, day == 'D2') %>%
mutate(sunrise = hms(sunrise),
sunset = hms(sunset))
#-- Prepare summary data for plotting --#
actogram_summary <- actograms %>%
# Select the required columns
select(year, month, data) %>%
# Unnest the data
unnest() %>%
# Group by year/month/day
group_by(year, month, day) %>%
# Calculate mean daily temperature
mutate(mean_temp = mean(body_temp)) %>%
# Calculate difference between 5-min temperatures and mean daily temperature
mutate(delta_temp = body_temp - mean_temp) %>%
# Re-group by year/month/time of day
group_by(year, month, time) %>%
# Calculate mean at each time point across each month for each year
summarise(summary_temp = mean(delta_temp)) %>%
ungroup() %>%
# Format year/month/time columns
mutate(month = fct_relevel(factor(month),
'January', 'February', 'March', 'April',
'May', 'June', 'July', 'August', 'September',
'October', 'November', 'December'),
time = hms(time))
#-- Plot --#
ggplot() +
geom_rect(data = photo_D1,
aes(xmin = sunrise,
xmax = sunset),
ymin = -Inf,
ymax = Inf,
fill = '#FFD700',
colour = '#FFD700',
size = 1,
alpha = 0.4) +
geom_rect(data = photo_D2,
aes(xmin = sunrise,
xmax = sunset),
ymin = -Inf,
ymax = Inf,
fill = '#A89070',
colour = '#A89070',
size = 1,
alpha = 0.4) +
geom_area(data = actogram_summary,
aes(x = time,
y = summary_temp,
fill = year),
colour = '#000000',
size = 0.2,
position = 'identity',
alpha = 0.6) +
scale_x_time(breaks = seq(from = 0,
to = 86400,
by = 21600), # Time measured in seconds
labels = c('00:00', '06:00', '12:00',
'18:00', '24:00')) +
scale_fill_manual(name = 'Year',
values = c('#FD7F28', '#2678B2', '#339F34')) +
scale_colour_manual(name = 'Year',
values = c('#FD7F28', '#2678B2', '#339F34')) +
labs(title = str_glue('{animal}: Average monthly difference in body temperature from\nmean daily temperature for each year'),
subtitle = 'Shaded area: Sunrise to sunset for the first (yellow) and last (tan) day of the month',
x = 'Time of day',
y = expression(Delta~body~temperature~from~daily~mean~(degree*C))) +
facet_grid(month ~ .,
drop = FALSE) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank(),
panel.spacing.y = unit(0, 'lines'))
Timing of daily minima and maxima
These analyses use hourly averages of 5-minute body temperature recordings. Median hourly temperature led to too many daily times with the same value and therefore multiple daily minima and maxima. So, we used the trimean, which produces a wider range of hourly temperature values, reducing the number of duplicate values for daily minima and maxima, while also avoiding skewing caused by acute excursions in body temperature that may occur if the arithmetic mean was used.
Density plots
Please note that the sequence of months is September 2017 (orange) through to March 2018 (blue).
#-- Define custom trimean function --#
trimean <- function(x, na = TRUE) {
Q25 <- quantile(x, na.rm = na, probs = 0.25)[[1]]
Q50 <- quantile(x, na.rm = na, probs = 0.50)[[1]]
Q75 <- quantile(x, na.rm = na, probs = 0.75)[[1]]
TM <- (2 * Q50) + Q25 + Q75
TM <- TM / 4
return(TM)
}
#-- Calculate hourly trimean temp (slow)--#
data_hour <- data %>%
# Round minutes to the hour (floor)
mutate(date_hour = floor_date(date_time, unit = 'hour')) %>%
# Group by data_hour and calculate trimean
group_by(date_hour) %>%
summarise(tri_mean = trimean(body_temp)) %>%
ungroup()
#-- Find time of the minimum/maximum hourly trimean temp for each day --#
# Minima
data_min <- data_hour %>%
separate(col = date_hour,
into = c('date', 'hour'),
sep = ' ',
remove = FALSE) %>%
group_by(date) %>%
filter(tri_mean == min(tri_mean,
na.rm = TRUE)) %>%
ungroup()
# Maxima
data_max <- data_hour %>%
separate(col = date_hour,
into = c('date', 'hour'),
sep = ' ',
remove = FALSE) %>%
group_by(date) %>%
filter(tri_mean == max(tri_mean,
na.rm = TRUE)) %>%
ungroup()
#-- Are there days with multiple minima/maxima? --#
# Minima
multi_min <- data_min %>%
group_by(date) %>%
# Get days with multiple minima
summarise(count = n()) %>%
filter(count > 1) %>%
ungroup() %>%
# Make them pretty
unite(col = 'days',
date, count,
sep = ' (n = ') %>%
mutate(days = paste0(days, ')')) %>%
.$days %>%
paste(., collapse = ', ') %>%
str_replace_all(pattern = '(.{39})\\s', # stuff in () is a group
replacement = '\\1\n') # \1 replaces the group
# Maxima
multi_max <- data_max %>%
group_by(date) %>%
# Get days with multiple minima
summarise(count = n()) %>%
filter(count > 1) %>%
ungroup() %>%
# Make them pretty
unite(col = 'days',
date, count,
sep = ' (n = ') %>%
mutate(days = paste0(days, ')')) %>%
.$days %>%
paste(., collapse = ', ') %>%
str_replace_all(pattern = '(.{39})\\s', # stuff in () is a group
replacement = '\\1\n') # \1 replaces the group
#-- Prepare data_min/max for plotting --#
# Minima
data_min %<>%
# Extract year, month and hour of day
mutate(year = year(date),
year = as.character(year)) %>%
mutate(month = month(date,
label = TRUE,
abbr = FALSE),
month = fct_rev(month)) %>%
mutate(hr = hour(date_hour))
# Maxima
data_max %<>%
# Extract year, month and hour of day
mutate(year = year(date),
year = as.character(year)) %>%
mutate(month = month(date,
label = TRUE,
abbr = FALSE),
month = fct_rev(month)) %>%
mutate(hr = hour(date_hour))
#-- Generate plot objects --#
# Minima
ggplot(data = data_min) +
aes(x = hr,
y = month,
fill = year,
colour = year) +
geom_density_ridges2(scale = 0.95,
alpha = 0.6) +
scale_x_continuous(breaks = c(0, 3, 6, 9, 12, 15, 18, 21),
labels = c('00:00', '03:00', '06:00',
'09:00', '12:00', '15:00',
'18:00', '21:00'),
limits = c(0, 23)) +
scale_fill_manual(name = 'Year',
values = c('#FD7F28', '#2678B2', '#339F34')) +
scale_colour_manual(name = 'Year',
values = c('#FD7F28', '#2678B2', '#339F34')) +
labs(title = str_glue('{animal}: Monthly density plot of time of body temperature minima'),
subtitle = '(Minima extracted from trimean hourly body temperature data)',
caption = str_glue('Days with multiple minima:\n {multi_min}'),
x = 'Hour of the day',
y = 'Month') +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank(),
plot.caption = element_text(size = 10,
hjust = 0))
# Maxima
ggplot(data = data_max) +
aes(x = hr,
y = month,
fill = year,
colour = year) +
geom_density_ridges2(scale = 0.95,
alpha = 0.6) +
scale_x_continuous(breaks = c(0, 3, 6, 9, 12, 15, 18, 21),
labels = c('00:00', '03:00', '06:00',
'09:00', '12:00', '15:00',
'18:00', '21:00'),
limits = c(0, 23)) +
scale_fill_manual(name = 'Year',
values = c('#FD7F28', '#2678B2', '#339F34')) +
scale_colour_manual(name = 'Year',
values = c('#FD7F28', '#2678B2', '#339F34')) +
labs(title = str_glue('{animal}: Monthly density plot of time of body temperature maxima'),
subtitle = '(Maxima extracted from trimean hourly body temperature data)',
caption = str_glue('Days with multiple minima:\n {multi_max}'),
x = 'Hour of the day',
y = 'Month') +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank(),
plot.caption = element_text(size = 10,
hjust = 0))
Heatmap
Unlike the density plots, the days with multiple minima or maxima were excluded from the heatmap analysis.
Please note that the sequence of months is September 2017 through to March 2018.
#-- Label maxima and minima data --#
# Minima
data_min %<>%
mutate(value = 'minimum')
# Maxima
data_max %<>%
mutate(value = 'maximum')
#-- Filter out days with multiple max/min values --#
# The spread function (used later) does not cope with duplicates
# Fix multi_max/min for use as a filter
## Minima
filter_min <- str_split(string = multi_min,
pattern = '\\s\\(.{5}\\)[,]\\s')[[1]] %>%
str_replace(pattern = '\\s\\(.{5}\\)',
replacement = '')
## Maxima
filter_max <- str_split(string = multi_max,
pattern = '\\s\\(.{5}\\)[,]\\s')[[1]] %>%
str_replace(pattern = '\\s\\(.{5}\\)',
replacement = '')
#-- Filter dataframes --#
# Minima
min_filtered <- data_min %>%
filter(!date %in% filter_min)
# Maxima
max_filtered <- data_max %>%
filter(!date %in% filter_max)
#-- Generate vector of removed days --#
multi_minmax <- c(filter_min, filter_max) %>%
# Filter for strings on length 10
.[str_detect(., '.{10}')] %>%
# Extract unique values
unique(.) %>%
# Sort by date
sort(.) %>%
# Collapse into a single string
paste(., collapse = ', ') %>%
# Add line-breaks
str_replace_all(pattern = '(.{39})\\s',
replacement = '\\1\n')
#-- Join min and max dataframes and organise data --#
data_minmax <- min_filtered %>%
bind_rows(max_filtered) %>%
# clean-up columns
select(date, hour, month, value) %>%
# Spread value column
spread(key = value,
value = hour) %>%
# Retain complete cases only
filter(complete.cases(.))
#-- Calculate time in minutes between min and max --#
data_minmax %<>%
mutate(interval = interval(paste(date, minimum),
paste(date, maximum)),
length_min = interval / dminutes(1))
#-- Plot --#
data_minmax %>%
select(-interval) %>%
mutate(maximum = hour(hms(maximum)),
minimum = hour(hms(minimum)),
month = fct_rev(month)) %>%
arrange(date) %>%
ggplot(data = .) +
aes(x = minimum,
y = maximum) +
geom_hex(bins = 12) +
geom_abline(intercept = 0,
slope = 1,
colour = '#999999') +
scale_y_continuous(limits = c(0, 24),
expand = c(0, 0),
breaks = seq(0, 24, 4)) +
scale_x_continuous(limits = c(0, 24),
expand = c(0, 0),
breaks = seq(0, 24, 4)) +
scale_fill_viridis_c(name = 'Days') +
labs(title = str_glue('{animal}: Hexagonal heatmap of time of daily minima and maxima'),
subtitle = 'Black diagonal line: line of identity',
caption = str_glue('Removed days with multiple maxima or minima: \n{multi_minmax}'),
x = 'Time of day of minimum (hours)',
y = 'Time of day of maximum (hours)') +
facet_wrap(~ month,
ncol = 3,
drop = FALSE) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank(),
plot.caption = element_text(size = 10,
hjust = 0))
Time between daily minimum and maximum
Unlike the density plots, the days with multiple minima or maxima were excluded from the heatmap analysis.
Please note that the sequence of months is September 2017 (orange) through to March 2018 (blue).
#-- Generate plot --#
# Interval calculated in previous code block
data_minmax %>%
# Extract year
mutate(year = year(date),
year = as.character(year)) %>%
ggplot(data = .) +
aes(x = length_min,
y = month,
colour = year,
fill = year) +
geom_density_ridges2(scale = 0.95,
alpha = 0.6) +
geom_vline(xintercept = 0,
colour = 'red',
size = 1) +
labs(title = str_glue('{animal}: Density plot of the time difference between\ndaily minimum and maximum body temprature'),
caption = str_glue('Removed days with multiple maxima or minima:\n {multi_minmax}'),
x = 'Time difference (min - max; minutes)',
y = 'Season') +
scale_fill_manual(name = 'Year',
values = c('#FD7F28', '#2678B2', '#339F34')) +
scale_colour_manual(name = 'Year',
values = c('#FD7F28', '#2678B2', '#339F34')) +
scale_x_continuous(limits = c(-1800, 1800),
breaks = seq(-1800, 1800, 600)) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank(),
plot.caption = element_text(size = 10,
hjust = 0))
Periodogram
Periodograms provide an estimate of the spectral density of a signal, with the assumption being that if there is a circadian rhythm, there will be a peak in density around a frequency of 24 hours. To limit the impact that changes in photoperiod across the year may have on the phase of the rhythm, we calculated monthly periodograms rather than assessing the complete data record. We used the multitaper method of signalling processing, with seven tapers and a time-bandwidth parameter of four.
#-- Nest data by year and month--#
periodgrams <- data %>%
mutate(month = month(date_time,
label = TRUE,
abbr = FALSE)) %>%
mutate(month = as.character(month)) %>%
group_by(year, month) %>%
nest() %>%
mutate(month = factor(month,
levels = c('January', 'February', 'March', 'April',
'May', 'June', 'July', 'August', 'September',
'October', 'November', 'December',
ordered = TRUE)))
#-- Generate periodogram data for each month --#
periodgrams %<>%
# Perform spectral analysis of the time series
# using the multitaper method (using default settings: nw = 4, k = 7)
mutate(periodogram = map(.x = data,
~ spec.mtm(ts(.x$body_temp),
nw = 4, k = 7,
dtUnits = 'day',
deltat = 1/288,
jackknife = TRUE,
plot = FALSE,
log = 'yes'))) %>%
# Add freq, spec, and jacknife 95% CIs into a dataframe
mutate(periodogram_df = map(.x = periodogram,
~ tibble(frequency = .x$freq,
spectrum = .x$spec,
upper_CI = .x$mtm$jk$upperCI,
lower_CI = .x$mtm$jk$lowerCI)))
#-- Generate periodgram and F-value plots for each month --#
## Days per month
days <- map(.x = periodgrams$data,
~ round(nrow(.x) / 288))
## Plot
periodgrams %<>%
# Periodogram plots
mutate(periodogram_plot = pmap(.l = list(periodogram_df,
as.character(month),
year, animal, days),
~ if(..5 < 16) {
ggplot(data = ..1) +
aes(x = frequency,
y = spectrum) +
geom_blank() +
labs(title = str_glue('{..4}: {..2} {..3} periodogram'),
subtitle = str_glue('No plot: series length < 16 days'),
x = 'Period (days)',
y = 'Spectral density') +
scale_x_continuous(limits = c(0, 5),
expand = c(0, 0),
breaks = c(0.01, 0.5, 1,
2, 3, 4, 5),
labels = as.character(round(c(1/0.01, 1/0.5, 1/1, 1/2, 1/3, 1/4, 1/5), 2))) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank())
} else {
ggplot(data = ..1) +
aes(x = frequency,
y = spectrum) +
geom_ribbon(aes(ymin = lower_CI,
ymax = upper_CI),
fill = 'wheat2',
colour = 'wheat3') +
geom_line() +
labs(title = str_glue('{..4}: {..2} {..3} periodogram'),
subtitle = str_glue('Series length (days): {..5} \nShaded area: 95% jackknife confidence interval \nMultitaper method: time-bandwidth (nw) = 4, tapers (k) = 7'),
x = 'Period (days)',
y = 'Spectral density') +
scale_x_continuous(limits = c(0, 5),
expand = c(0, 0),
breaks = c(0.01, 0.5, 1,
2, 3, 4, 5),
labels = as.character(round(c(1/0.01, 1/0.5, 1/1, 1/2, 1/3, 1/4, 1/5), 2))) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank())}))
# Print plots
walk(periodgrams$periodogram_plot,
~print(.x))
Autocorrelation
Like the periodogram analysis, to limit the impact that changes in photoperiod across the year may have on the phase of the rhythm, we calculated autocorrelation for monthly data rather than assessing the complete data record.
In addition, we reduced the size and smoothed the data before performing the autocorrelation by averaging the 5-minute body temperature data across 48 30-minute duration bins per day. The autocorrelation was performed with a maximum lag of 480 30-minute bins (10 days).
#-- Nest data by year and month--#
autocorr <- data %>%
mutate(month = month(date_time,
label = TRUE,
abbr = FALSE)) %>%
mutate(month = as.character(month)) %>%
group_by(year, month) %>%
nest() %>%
mutate(month = factor(month,
levels = c('January', 'February', 'March', 'April',
'May', 'June', 'July', 'August', 'September',
'October', 'November', 'December',
ordered = TRUE)))
# Process data
autocorr %<>%
# Prepare for plotting
mutate(month = as.character(month)) %>%
# Create 30-minute binned body temperature data
mutate(temp_binned = map(.x = data,
~ .x %>%
# Extract minutes
mutate(time_min = minute(date_time),
time_hour = hour(date_time)) %>%
# Re-level minutes to 30 minute intervals
mutate(time_30 = ifelse(time_min > 0 & time_min < 35,
yes = '00:00',
no = '30:00'),
time_30 = str_glue('{time_hour}:{time_30}')) %>%
# Bin by 30-min interval and calculate
# average temperature per bin
group_by(day, time_30) %>%
summarise(temp_binned = mean(body_temp, na.rm = TRUE)) %>%
ungroup() %>%
# Extract body_binned
.$temp_binned)) %>%
# Generate acf object (autocorrelation lag.max = 10 days (n = 480))
mutate(ACF = map(.x = temp_binned,
~ acf(.x,
lag.max = 480,
plot = FALSE))) %>%
# Is the length of temp_binned less than 10 days (480 bins)
mutate(data_length = map(.x = temp_binned,
~ ifelse(length(.x) < 480,
yes = 'TRUE',
no = 'FALSE'))) %>%
# Plot data
mutate(acf_plots = pmap(.l = list(ACF, as.character(month), year, animal, data_length),
~ if(..5 == 'FALSE') {
autoplot(..1,
conf.int.fill = '#2678B2') +
labs(title = str_glue('{..4}: {..2} {..3} autocorrelation of body temperature'),
subtitle = 'Temperature data averaged across 30-minute bins \nMaximum lag: 10 days (480 bins) \nBlue ribbon: 95% confidence interval of a white noise process',
x = 'Days',
y = 'Autocorrelation coefficient') +
scale_y_continuous(breaks = seq(from = -0.4,
to = 1,
by = 0.2),
limits = c(-0.4, 1.1),
expand = c(0, 0)) +
scale_x_continuous(breaks = seq(from = 0,
to = 480,
by = 48),
labels = 0:10,
limits = c(-5, 485),
expand = c(0,0)) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank())
} else {
df <- data.frame(acf = ..1$acf[1:nrow(..1$acf), , 1],
lag = ..1$lag[1:nrow(..1$lag), , 1])
ggplot(data = df) +
aes(x = lag,
y = acf) +
geom_blank() +
labs(title = str_glue('{..4}: {..2} {..3} autocorrelation of body temperature'),
subtitle = 'No Plot: Number of days < maximum lag period (10 days)',
x = 'Days',
y = 'Autocorrelation coefficient') +
scale_y_continuous(breaks = seq(from = -0.4,
to = 1,
by = 0.2),
limits = c(-0.4, 1.1),
expand = c(0, 0)) +
scale_x_continuous(breaks = seq(from = 0,
to = 480,
by = 48),
labels = 0:10,
limits = c(-5, 485),
expand = c(0,0)) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank())}))
# Print plots
walk(autocorr$acf_plots,
~print(.x))
Session information
sessionInfo()## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] hexbin_1.27.2 zoo_1.8-4 skimr_1.0.5
## [4] multitaper_1.0-14 lubridate_1.7.4 ggridges_0.5.1
## [7] ggfortify_0.4.5 forcats_0.4.0 stringr_1.4.0
## [10] dplyr_0.8.0.1 purrr_0.3.1 readr_1.3.1
## [13] tidyr_0.8.3 tibble_2.1.1 ggplot2_3.1.0
## [16] tidyverse_1.2.1 magrittr_1.5
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 xfun_0.5 reshape2_1.4.3
## [4] haven_2.1.0 lattice_0.20-38 colorspace_1.4-1
## [7] generics_0.0.2 viridisLite_0.3.0 htmltools_0.3.6
## [10] yaml_2.2.0 rlang_0.3.2 pillar_1.3.1
## [13] glue_1.3.1 withr_2.1.2.9000 modelr_0.1.4
## [16] readxl_1.3.0 plyr_1.8.4 munsell_0.5.0
## [19] gtable_0.3.0 cellranger_1.1.0 rvest_0.3.2
## [22] evaluate_0.13 labeling_0.3 knitr_1.21
## [25] highr_0.7 broom_0.5.1 Rcpp_1.0.1
## [28] scales_1.0.0 backports_1.1.3 jsonlite_1.6
## [31] gridExtra_2.3 hms_0.4.2 digest_0.6.18
## [34] stringi_1.4.3 grid_3.5.2 cli_1.1.0
## [37] tools_3.5.2 lazyeval_0.2.2 crayon_1.3.4
## [40] pkgconfig_2.0.2 xml2_1.2.0 assertthat_0.2.1
## [43] rmarkdown_1.11 httr_1.4.0 rstudioapi_0.9.0
## [46] R6_2.4.0 nlme_3.1-137 compiler_3.5.2